Time-dependent transport in graphene nanoribbons 
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We theoretically investigate the time-dependent ballistic transport in metallic graphene nanorib- 
bons after the sudden switch-on of a bias voltage V. The ribbon is divided in three different regions, 
namely two semi-infinite graphenic leads and a central part of length L, across which the bias drops 
linearly and where the current is calculated. We show that during the early transient time the 
system behaves like a graphene bulk under the influence of a uniform electric field E — V/L. In 
the undoped system the current does not grow linearly in time but remarkably reaches a temporary 
plateau with dc conductivity <7i = ire' 2 /2h, which coincides with the minimal conductivity of two- 
dimensional graphene. After a time of order L/vf (vf being the Fermi velocity) the current departs 
from the first plateau and saturates at its final steady state value with conductivity 02 = 2e 2 /h 
typical of metallic nanoribbons of finite width. 

PACS numbers: 



The recent isolation of single layers of carbon atoms^ 
has attracted growing attention in the transport prop- 
erties of graphene-based devices. In these systems un- 
conventional phenomena like the half-integer quantum 
Hall effectPl and the Klein tunneling^ have been ob- 
served. Such peculiar behavior stems from the relativistic 
character of the electrons in the carbon honeycomb lat- 
tice. Closed to Dirac point the charge carriers behave 
as two-dimensional (2D) massless Dirac fermions^ and 
have very high mobility^. This fact has stimulated the- 
oretical and experimental investigations into graphenic 
ultrafast devices^ like field effect transistors^, p-n junc- 
tion diodes and THz detectors^]. In these systems it 
is crucial to have full control of the electronic response 
after the sudden switch-on of an external perturbation, 
and the study of the real-time dynamics is becoming 
increasingly important. The investigation of the tran- 
sient response also is of fundamental interest. One of 
the most debated aspects of the transport properties of 
graphene is the minimum conductivity cr m j n at the Dirac 
point. From the theoretical point of view the problem 
arises from the fact that the value of <7 m i n is sensitive 
to the order in which certain limits (zero disorder and 
zero frequency) are takerP^, th us produci ng different val- 
ues around the quantum e 2 /^ 10 * 11 * 12 * 13 !. Very recently 
Lewkowicz and Rosenten overcame this ambiguity em- 
ploying a time-dependent approach in which the dc con- 
ductivity is calculated by solving the quench dynamics of 
2D Dirac excitations after the sudden switching of a con- 
stant electric field. Interestingly their approach does not 
suffer from the use of any regularization related to the 
Kubo or Landauer formalism and yields <r min = 7re 2 /2h. 
Despite the large effort devoted to the study of the trans- 
port properties of graphenic systems, a genuine real-time 
analysis which treats on equal footing transients effects 
and the long-time response of graphene nanoribbons in 
contact with semi-infinite reservoirs is still missing. 



In this Letter we study the time-dependent trans- 
port properties of undoped graphene nanoribbons with 
finite width and virtually infinite length after the sud- 
den switch-on of an external bias voltage. The geometry 
that we consider is sketched in Fig[T] The nanoribbon 
is divided in three regions, namely a left (L) and a right 
(R) semi-infinite graphenic reservoirs and a central (C) 
region of length L = aN c , where N c is the number of 
cells along the longitudinal x direction and a = 2.46 A is 
the graphene lattice constant. The width of the ribbon is 
W = a\/3Ny, where N y is the number of cells along the 
transverse y directi on, in w hich periodic boundary con- 
ditions are imposed 15 * 16 * 17 !. The three regions are linked 
via transparent interfaces, in such a way that in equilib- 
rium the system is translationally invariant along the x 
direction, see Fig[T] 

Once the system is driven out of equilibrium, the 
quench dynamics involves high energy excitations and 
the Dirac-cone approximation, which is valid only at en- 
ergies lower than 1 eV, is inaccurate. For this reason 
we adopt a tight-binding description of the system, with 
Hamiltonian given by 

H = H a + U{t) = v J2 4 c i + 6(t) VicU , (!) 

where the spin index has been omitted and v = 2.7 eV is 
the hopping integral of graphene. The first sum runs over 
all the pairs of nearest neighbor sites of the ribbon hon- 
eycomb lattice and is the annihilation (creation) op- 
erator of a 7r electron on site i. Here we use the collective 
index i = {p, i yi i x } to identify a site in the nanoribbon 
such that p — a, (3 indicates the two inequivalent longitu- 
dinal zig-zag chains, i y denotes the cell in the y direction, 
and i x is the position in the x direction. Hq describes the 
translationally invariant equilibrium system, while U(t) 
is the bias perturbation with (non self-consistent) voltage 
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FIG. 1: Schematic representation of the system. The nanorib- 
bon has semi-infinite L a nd R graphenic reservoirs and peri- 
odic boundary conditions are imposed along the y direction. 
The bias voltage profile with a linear drop inside the C region 
is also shown. 
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where V is the total applied voltage and i x G (0, L) is the 
x-coordinate of the site i in the C region, which is subject 
to the uniform electric field E = V/L (see Fig(T]). The 
above modelling of the bias profile could find an approx- 
imate realization, e.g., in a planar junction in which the 
L and R regions of the nanoribbon are on top of metallic 
electrodes. 

The time-dependent total current I(t) flowing across 
the interface in the middle of the C region is written as 



*(*) = 2E E J ? } W' 



(3) 



i y =l p=a,/3 



where the factor 2 accounts for the spin degeneracy and 
ij^ is the current flowing across the p = a, (3 chain in 
the iy-th cell, in the middle of the C region (see Figll]) . 
Since periodic boundary conditions are imposed along 

(%>) 

the y direction, the current J,- does not depend on the 
cell and chain indices i y and p, and hence I(t) — AN y I(t), 
where we have defined I = J„- = /„• for any i v . Since 

Ly ly •■> y 

the transverse momentum k y = 27m/ 'y/3aN y (with n 



0, 



. N y — 1) is conserved the current / can be written as 



I(t) 



(4) 



It is seen that the total current of the original N y - 
wide problem can be calculated by summing the cur- 
rents Jfe coming from N y independent 1-wide ladder 
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FIG. 2: Ladder model for fixed transverse momentum k y . 
The bias profile is the same as in Figjl] 

problems^ with staggered fc^-dependent transverse hop- 
ping (see Fig{2|, with the same bias profile as in Eq.Q. 

In order to calculate Ik it is convenient to introduce 
the y-Fourier transform of the original electron operators 
C{^ H ,U = (^)"' /2 E iy e i& v%^ C{pA ^ } in terms of 
which the current 1^ can be cast as 
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where G < is the lesser Keldysh Green's function 
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The time-evolution of the Green's function is evaluated 
according to 
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I {p,k v ,i x };{p,k v ,j x } 



(7) 



where we recall that there is no actual dependence on p 
and where / is the Fermi distribution function. For each 
transverse momentum k y the current Ik is numerically 
calculated by computing the exact time evolution of the 
corresponding ladder system in Fig(2j where we take the 
reservoirs with a finite length L r — aN r . This approach 
allows us to reproduce the time evolution of the infinite- 
leads system up to a time T max w 2L r /vp, where vp is 
the Fermi velocitj^. For t > T max electrons have time 
to propagate till the far boundary of the leads and back, 
yielding undesired finite size effects in the calculated cur- 
rent. Accordingly we choose N r such that T max is much 
larger than the time at which the steady-state is reached. 

For practical purposes we represent Hq and H{t) in a 
hybrid basis, in which the Hamiltonians of the isolated 
(equilibrium and biased) L and R reservoirs are diag- 
onal in the set of states derived analytically in Ref.22 
while the Hamiltonian describing the C region remains 
represented in the basis {p, k y ,i x }. The enormous advan- 
tage of such choice is clarified in the following. Since we 
are interested in the dc conductance we set a small bias 
V « hvF/eW . In this way we ensure that at long times 
a linear regime is established in which only the electrons 
right at the Dirac point (with k y = 0) contribute to the 
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FIG. 3: Total time dependent current I(t) with geometric 
parameters L ~ 25 nm, W ~ 106 nm and applied volt- 
age V = 8 x 1CP 4 Volt. Current is in units of Voo, where 
(Jo = e 2 //i is the quantum of conductance and time in units 
of t v — 2.5 x 10~ 16 s. The dashed line represents the 
long-time conductance 02/0-0 = 4 given by the Landauer 
fbrmulsP, while the dotted line represents the combination 

(Omin X W/L)/O = 7r/2 X W / L. 



total current. This agrees with the Landauer formula, 
according to which only the states within the bias win- 
dow contribute at the steady-state. On the other hand 
during the transient all transverse modes are excited and 
the sum in Eq.Q must be computed including the com- 
plete set of k y . Indeed we have checked numerically that 
in the limit t — > 00 all currents Jt , except the one with 
k y = 0, vanish. We have also observed that the damping 
time of Ik y >o(t) goes like ky 1 and the calculation of the 
currents with small k y requires a very long propagation 
before the zero-current steady-state is approached. As a 
consequence very large values of N r are needed, making 
the computation in principle too demanding. Such nu- 
merical difficulty is, however, compensated by the fact 
that, after a transient time of order L/vp, for k y close to 
the Dirac point only few low-energy states contribute to 
Jfe (t). Therefore for any given k y > we introduce an 
energy cutoff Ak y ~ I0hvpk y in the reservoirs Hamiltoni- 
ans and retain only the lead-eigenstates with transverese 
momentum k y and energy in the range (— Ak y , At). This 
allows us to deal with very long leads (N r >> 1) and can 
be explicitly implemented by using the analytic eigen- 
states of a rectangular grapheme macromolecule^l. The 
C region (which is the one where we calculate the current) 
is treated exactly within the original full basis {p, k y ,i x }, 
but the overall computational cost remains moderate. 

In Fig{3]we show the time-dependent current I(t) cal- 
culated for a nanoribbon of central length L ss 25 nm 
(N c = 100) and width W w 106 nm (N y = 250) with an 
applied voltage V — 8 x 10^ 4 Volt and zero temperature. 
The very early transient regime (t < h/v = t v ) depends 
on the details with which the electric field E has been 
switched-on (sudden in time in the present case). For 
ty ^ t < L/vf, however the time evolution only depends 
on the geometry of the device and on the potential profile. 
Within this time domain, if the ribbon is wide enough, 
the ballistic electron dynamics of the system probes only 



FIG. 4: Short time transient current I(t) for different ribbon 
widths W\ = 106 nm, W 2 = 27 nm, W 3 = 13 nm. The rest 
of parameters, dotted/dashed lines and units are as in Fig[3] 

the bulk properties of graphene, since the particles do 
not have time to explore the reservoirs, where E = 0. 
According to the Drude picture, the ballistic transport 
in bulk materials subjected to uniform electric fields pro- 
duces a time-dependent conductance given by 



where 7 is a constant proportional to the density of 
states. In ordinary solids the dc conductance (ui — > 0) 
increases linearly in time up to the breakdown of the 
ballistic regime, in which one has to replace t with the 
finite scattering time r. However in pure graphene the 
transport is ballistic up to very large length-scales of the 
order of micron and a can apparently diverge, produc- 
ing a Drude peak. Nevertheless in undoped graphenic 
samples the density of states at the Fermi level vanishes, 
thus making the product yt constant at long times. This 
subtle compensation is at the origin of the finite minimal 
dc conductivity of 2D pure graphene^ and of the diffi- 
culties in constructing a suitable propagation scheme in 
finite width nanoribbons. In Fig(3]we see that W w 100 
nm is enough to observe such compensation. The cur- 
rent I(t), instead of increasing linearly in time, reaches 
a temporary plateau with average current I\ which lasts 
till t sa L/vp ss 60t„. On the contrary in an ordinary 
system (e.g. in a 2D square lattice model) the current 
would be linear in time up to L/vp, producing the well 
known Drude peak in the bulk limit L, W —> 00 at fixed 
E. We would like to observe that the plateau at L is 
reached via transient oscillations with frequency 2v/$rQ 
Interestingly such frequency is not displayed by any of the 
individual currents y (t), but appears only as a cumu- 
lative effect after the summation in Eq.Q is performed. 
Therefore it is a a genuine bulk property and may be 
at the origin of the resonant effect predicted to occur in 
optical response of graphene right at u> = 2v/tfr^. 

From the first plateau of I(t) we can provide an in- 
dependent evaluation of the minimal conductivity of 
graphene. According to its definition, the conductivity 
(Ti of a bulk system subjected to a small constant elec- 
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FIG. 5: Short time transient current I(t) for three different 
metallic armchair ribbons with open boundary conditions. 
The ribbon parameters are L — 5.5 nm, Wi = 3.6 nm, 
W2 = 2.1 nm, Wz = 0.6 nm and the applied bias is V = 0.03 
Volt. Units and dotted lines are as in Fig[3] The dashed line 
represents the lo ng-tim e conductance a% /(To = 2 given by the 
Landauer formula^ 11 * 24 !. 



trie field E is given by the current density J divided by 
E: 



°i = ^ = 



h L 
V W 



(9) 



By exploiting Eq.([9]) it is seen that our data are consis- 
tent with the value o\ — ire 2 /2h = <7 m i n with excellent 
precision, see Fig|3] 

Further we have investigated the formation of the first 
plateau as a function of the ribbon width W. In Figj4] 
we see that for narrow ribbons with W ~ 10 nm finite- 
size effects produce a drastic deviation from the ideal 
graphene bulk. The transient current does not show a 
temporary saturation, but grows with an approximate 
linear envelope (up to a time ~ L/vf), in qualitative 
agreement with the Drude behavior given in Eq.([8|. 

At times larger than ~ L/vp the electrons start explor- 
ing the reservoirs, where the electric field is zero. At this 
point a standard dephasing mechanism sets in and the 
current tends to its true final steady-state. As discussed 
above, however, such value is reached after a very slow 
damping process, in which the current displays decaying 
oscillations with dominant frequency Q = 2nvp/W (see 
Figj3|, hw being the energy spacing between the trans- 
verse energy subbands of the ribbon. We have checked 
numerically that the asymptotic value / 2 agrees well with 
the Landauer formula and does not depend on L and W , 



provided that V << Hvp/eW. Thus the conductance of 
the device is simply extracted as &2 = I2/V. Our nu- 
merical data provides 172 = 4e 2 //i with high numerical 
accuracy. This value is indeed twice the conductance of 
metallic nanoribbons, and this is due to the chosen pe- 
riodic boundary conditions^. Thus we have calculated 
I{t) also in the case of open boundary conditions. Un- 
fortunately, as discussed above, the lack of translational 
invariance along the transverse direction makes the com- 
putation much more demanding, and only system with 
small L and W can be studied within the present ap- 
proach. In Fig{5 we show I{t) for three different metallic 
armchair nanoribbons with open boundaries. It can be 
seen that already for W ~ 2 — 4 nm there is a tendency 
to form the universal first plateau leading to er m ; n , while 
for W < 1 nm the current I(t) grows linearly in time 
until t ss L/vp. On the other hand at long times the 
current tends clearly to the Landauer value, cons istent 
with 02 = 2e 2 //i, independently on the aspect raticS32U. 

In summary we pointed out the subtle difficulties 
in constructing a reliable method to perform time- 
evolutions of finite width graphene nanoribbons and pro- 
posed an efficient numerical scheme to overcome them. 
We presented a real-time study of the transport proper- 
ties of these systems in contact with virtual semi- infinite 
reservoirs in the linear regime. We have shown that 
for large enough undoped samples the time-dependent 
current displays two plateaus. From the first of these 
plateaus we can extract an independent measure of the 
minimal conductivity ire 2 /2 of bulk graphene by resort- 
ing the aspect ratio L/W of the device. The second 
plateau corresponds to reaching the steady-state and 
is independent of the geometry. Here the conductance 
is 2e 2 /h, which coincides with the Landauer result for 
metallic nanoribbons. To conclude we wish to point 
out that in presence of ac bias, the time-dependent con- 
ductivity can be used to obtain the optical conductivity 
<j ac of graphene. It was shown experimentally^ and ex- 
plained theoretically^ that cr ac is almost (^-independent 
and equals cr min with high accuracy over a wide range of 
frequencies. Remarkably to extract the universal value 
of <r ac high frequency signals with to ~ l/t v have been 
employecpSl Therefore we believe that a real-time ap- 
proach like to one presented here is needed to enlighten 
the crossover from the dc case to ultrafast scenarios in 
which the period of the ac signal is comparable with the 
intrinsic hopping time of the bulk system. 
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